Wavefunction correction scheme for non fixed-node diffusion Monte Carlo 
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Wavefunction correction scheme, which was developed as a variance reduction tool for the pure 
and fixed- node diffusion Monte Carlo (DMC) computations by Anderson and Freihaut, is applied to 
the DMC computations of fermions without using the fixed-node constraint. This technique is found 
to be suitable for the non fixed-node calculations because of the significant decreases observed in 
the computation times for calculating results with certain statistical error values in the benchmark 
computations. 



I. INTRODUCTION 



Antisymmetry condition of the wavefunction on iden- 
tical particle exchanges complicates the electronic struc- 
ture calculations of the fermionic many-body systems. 
The projector quantum Monte Carlo (QMC) methods 
[3, [IJ, facilitating imaginary time evolution of an ini- 
tial quantum state, are accurate tools for such fermionic 
calculations. However, the evolution in imaginary time 
tends to the symmetric bosonic ground state instead of 
the antisymmetric fermionic ground state, resulting in 
the famous fermion sign problem Most of the at- 
tempts for an exact imposition of the antisymmetry con- 
dition in the projector methods facilitate plus and mi- 
nus signed walkers diffusing and canceling each other 
whenever encounters occur d]. Usual population con- 
trol mechanism in these methods can only stabilize the 
difference between the plus and minus signed populations 
and therefore cancellations of opposite signed walkers are 
essential for controlling the total population. Population 
control problem arises for larger systems since the plus 
and minus signed walker encounter rate is very low in 
the higher dimensional configuration spaces. The cor- 
related dynamics of plus and minus signed walkers and 
the antisymmetric guiding functions used in the recently 
developed fermion Monte Carlo method [|| increase the 
cancellation rate to some extent. However, it was shown 
that these developments were not enough for the resolu- 
tion of the sign problem 6]. Therefore, applications of 
the exact methods are currently limited to a small num- 
ber of fermions. 

Wavefunction correction scheme was developed for the 
projector QMC computations as a variance reduction 
tool 0| . The difference between the true ground state 
wavefunction and a trial wavefunction is sampled in this 
technique instead of the ground state wavefunction itself. 
This technique is used as an efficiency improvement in the 
projector QMC computations of bosons [7-9] as well as 
fixed- node projector QMC computations of fermions [t|- 
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In the current study, wavefunction correction scheme 
is applied to plus-minus cancellation facilitating compu- 
tations of the fermionic systems without using the fixed- 
node approximation. Efficiency improvements due to the 
wavefunction correction scheme in such QMC calcula- 
tions are investigated on some simple benchmark sys- 
tems. Calculations are carried out using the diffusion 
Monte Carlo (DMC) [3, S, ESl method but the correction 
scheme is applicable to other projector QMC methods as 
well. 



II. METHOD OF COMPUTATION 

DMC, an highly accurate QMC method, relies on the 
fact that the form of the imaginary time Schrodinger 
equation is a diffusion equation with a source term: 



9 T *(x, r) = - V 2 *(x, t) - [V( X ) - £ fl ]*(x, r) , (1) 

where x is the position vector in the configuration space 
of the physical system. The potential energy V(x) defines 
the interactions between the particles and with external 
sources. DMC treats the wavefunction ^(x, r) as a den- 
sity distribution of some number of hypothetical parti- 
cles, also called as walkers, diffusing in the D x N dimen- 
sional configuration space, D being the number of space 
dimensions and N being the number of identical parti- 
cles. These walkers are subjected to a branching process 
according to the source term of the diffusion equation 
which is the term including the potential energy V(x) in 
the Schrodinger equation. Population control is estab- 
lished by controlling the rate of the branching process 
via adjustments of the reference energy En which is an 
overall energy shift. The evolution in imaginary time r 
projects out the ground state component of an arbitrary 
initial wavefunction ^(x, 0) in the long r limit [lO]. The 
DMC method uses short time propagator and thus the 
calculated expectation value result has a time step error 
which can be made insignificant by a time step extrapo- 
lation. 

When the fixed-node constraint is not enforced the sign 
problem manifests itself as the problem of imposing the 
antisymmetry condition. Minus signed walkers arise in 
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such non fixed-node calculations, breaking down the pop- 
ulation control mechanism for large systems even with 
the cancellations of opposite signed walkers. 

In the correction scheme, DMC method is modified 
to make correction on a known trial wavefunction 043- 
Wavefunction is divided into two parts as the trial wave- 
function $t(x) and the remaining unknown part $(x, r) 
which is sampled through the DMC calculation. It is 
helpful to substitute $(x, r) + $t(x) for the wavefunc- 
tion ^(x, t) in the imaginary time Schrodinger equation 
for comprehending the effect of the correction scheme: 

9 T $(x,r) = ±V 2 [$(x,r) + $ T (x)] 

(2) 

-[V(x)-.Er][$(x > t)+*t(x)] , 

which can be simplified using the definition of the local 
energy £x(x) = -H"$t(x)/$t(x) as follows: 

c> r $(x, t) = ±V 2 $(x, r) - [V{x) - £r]$(x, t) 

(3) 

-[E L {x)-E R ] $ t (x). 

Last term on the rhs of the above equation is the ad- 
ditional term related to the correction scheme whose 
sole effect may be simulated by some number of vac- 
uum branchings carried out along the simulation region 
with the branching factors linearly proportional to the 
extra term [£x(x) — Er] $t(x). A plus or minus signed 
walker depending on the sign of this term may be added 
to the walker population as a result of a single vacuum 
branching. In the current computations, some number 
of points in the configuration space are generated using 
the Metropolis algorithm according to the distribution 
3>t(x) which is positive definite in the simulation region 
described below and the branching factors are calculated 
using the factors [El (x) — Er] . 

Antisymmetry condition is imposed on the wavefunc- 
tion by using the concept of the permutation cell which 
is the repeating unit cell of the wavefunction for iden- 
tical fermions and bosons [ll|, Ell • Computation is car- 
ried out in a single permutation cell which is taken as 
a positive valued nodal region of the trial wavefunction 
$t( x ) (Trial wavefunctions are chosen to have nodal re- 
gions having permutation cell property. Density func- 
tional theory results also have this property [13]]). All 



the walkers of the plus and minus sign are initially gen- 
erated in the chosen permutation cell and the outgoing 
walkers are permuted back to the simulation region. A 
sign change is applied if an odd number of particle per- 
mutations are required to take the walker back inside the 
permutation cell. 

The normalization of the wavefunction is satisfied by 
keeping the plus and minus signed walker amounts equal 
to each other via step by step adjustments of the DMC 
reference energy which affects the rates of the walker 
and vacuum branchings. The trial wavefunction used for 
the correction is normalized separately to an optimum 
value. The ratio of the trial wavefunction normaliza- 
tion (J $T(x)<if2 integral is over the simulation region 
where $t(x) is positive) to the number of walkers from 
each sign is an important parameter (r n ) of the correc- 
tion scheme calculations. A larger value of this ratio in- 
creases the efficiency of the correction scheme calculation 
since the contribution of the walker distribution and the 
variance of the computed result decreases in such a case. 
However, if the ratio is adjusted to have a very high value 
the effect of the walker distribution gets insignificant and 
the wavefunction is not corrected. This mentioned ratio 
and consequently the efficiency of the method can be in- 
creased as the trial wavefunction gets closer to the true 
fermionic wavefunction. 

Population control is established by the cancellations 
of opposite signed walkers which encounter during the 
correlated random walk process in which the Gaussian 
random walk vectors of the paired walkers are symmet- 
ric with respect to the perpendicular bisector of the line 
connecting the two walkers Q. Cancellation process is 
temporarily stopped when the number of walkers de- 
creases below a certain threshold value for the small di- 
mensional calculations (up to four dimensional configu- 
ration spaces) since the cancellation rate is very high for 
them. 

Computation of the ground state energy expectation 
value should be modified accordingly. Necessary modifi- 
cations are derived by integrating both sides of the eigen- 
value equation over the simulation region volume. The 
final energy expression is: 



-\ J dn VgKx, T ).dS + /„ y(x)$(x, r)dSl + / n E L (x)* T (x)dSl 
r n [$(x,r) + $ T (x)]^ 



where divergence theorem is used to aquire the first term boundaries of the simulation region |12j | . Second term of 
of the numerator which is about the walker flow at the the numerator is the sum of walker potential energies in 
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the energy calculation of the usual DMC. Third integral 
of the numerator is about the trial wavefunction being 
corrected throughout the computation and it can be cal- 
culated in the beginning of the simulation using Monte 
Carlo integration technique without respecting the nor- 
malization of the trial wavefunction. The value of the 
integral f $T"(x)df2 of the denominator is given as a pa- 
rameter of the method (determines the parameter r n to- 
gether with the number of walkers from each sign) and 
the mentioned Monte Carlo integration is multiplied by 
this given value of the $t( x ) integral for the normaliza- 
tion issue. The remaining first integral of the denomi- 
nator is the difference between the number of plus and 
minus signed walkers. Ground state energy is calculated 
considering these separate terms in each time step and 
it is time averaged after some thermalization steps. The 
time step errors are ignored in the following benchmark 
calculations since they are very small compared to the 
statistical error bars. 



III. BENCHMARK COMPUTATIONS 

A. Harmonic fermions 

Harmonic fermions are preferred in the benchmark cal- 
culations for their property of being exactly solvable. An- 
alytical solutions are disturbed slightly to prepare trial 
wavefunctions suitable for testing the current method. 
Two fermion systems are studied which have Hamilto- 
nian functions in the following form: 



values are set to the highest values allowing the complete 
correction for each case since the efficiency of the correc- 
tion scheme increases as the r n increases. These values 
are easily determined since a deviation in the calculated 
energy expectation value starts to occur beyond a cer- 
tain point as the mentioned ratio increases. Calculated 
energy expectation value matches with the true value be- 
low this point which is identified as the optimum value 
of the parameter r n . 

Computations are also carried out without using the 
correction scheme for a comparison of the computational 
efforts of the two cases. DMC without any corrected 
trial wavefunctions is used for these comparison calcula- 
tions. Same permutation cells that used in the correction 
scheme computations are used where outgoing walkers 
are treated in the same way. Correlated walk of oppo- 
site signed walkers with the cancellation process is also 
facilitated in the comparison case computations. Com- 
putation times of the two cases for calculating the results 
with certain statistical error values are compared. 

Computation results for the harmonic fermion systems 
up to four space dimensions are given in TABLE U for 
the correction scheme computations. Energy expecta- 
tion values for the used trial wavefunctions (&r), cal- 
culated by Monte Carlo integration technique, are also 
given. Computed energies are very close to the true val- 
ues (Eqs) for the all cases. Efficiency improvements can 
be seen from the ratios of the comparison case computa- 
tion times to the correction scheme computation times. 
Significant decreases in the computation times are ob- 
served for the all studied space dimensions when the cor- 
rection scheme is used. 



H = + (rf+rl), (5) 

where the vectors ri and v 2 denote individual particle 
coordinates and lo 2 is a constant whose numerical value 
is taken as 0.03 in the current calculations. The trial 
wavefunction is chosen as: 

^ T = e ei ^ (r?+r ^(x2 +e 2 yl-xi- e 2 y 2 ), (6) 

where x, y are the particle coordinate components in two 
separate space dimensions and e±, e 2 are free parameters. 
This trial wavefunction gives the true fermionic ground 
state in £\ — > 1, e 2 — V limit for arbitrary number of 
space dimensions. A non zero value of the parameter 
e 2 distorts the nodal surface of the trial wavefunction 
for space dimensions larger than one. However, the per- 
mutation cell property of the nodal region is preserved 
which is taken as the simulation region where outgoing 
walkers arc taken inside as described in the previous sec- 
tion. The DMC time step is chosen as 0.003 dimension- 
less time units and data is collected for 80000 time steps 
after the thermalization steps. Number of points for the 
vacuum branching process is chosen as 500 and kept con- 
stant during the computations. The ratio parameter r n 



TABLE I: Correction scheme computation results for two har- 
monic fermions. d: space dimension, £i,£2: disturbance pa- 
rameter values, E c : calculated energy expectation value us- 
ing the correction scheme, Eq§: true value of the fermionic 
ground state energy, E^: trial wavefunction energy (all ener- 
gies are given in dimensionless units), N w : stabilized number 
of walkers from each sign, r n : ratio of the trial wavefunction 
normalization to the number of walkers from each sign, rt: ra- 
tio of the comparison case computation time to the correction 
scheme computation time. 



d £1 £2 E c E GS E T JVw r n n 

1 0.964 0.00 0.34644(31) 0.34641 0.34677 195 14.4 2.7 

2 1.000 0.05 0.51978(62) 0.51962 0.52205 495 5.7 5.5 

3 1.000 0.05 0.6932(11) 0.69282 0.70143 520 4.6 8.8 

4 1.000 0.05 0.8660(10) 0.86603 0.87473 1556 1.9 4.8 



Images for the trial wavefunction (top image) and its 
difference from the true fermionic wavefunction (middle 
image) is given in the FIG. Q] for the ID computation 
whose configuration space is two dimensional. Average 
walker distribution during the correction scheme DMC 
computation (bottom image) is also given. Walker dis- 
tribution is calculated in a single permutation cell and it 
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is reflected to the other cell with a sign inversion in order 
to generate a plot for the all configuration space. Minus 
signed walkers give negative weights when the average is 
calculated. Walker distribution fits well with the differ- 
ence function as expected when the correction scheme is 
used. 



6 r 




FIG. 1: Wavefunction images for correction scheme DMC 
computation of the two harmonic fermions in ID. TOP: Trial 
wavefunction. MIDDLE: Difference between the trial wave- 
function and the true fermionic ground state. BOTTOM: 
Average walker distribution during the DMC computation. 



B. Helium atom lowest triplet state: ls2s 3 S 

As a more realistic example, Helium atom lowest 
triplet state energy eigenvalue is computed using the 
current method. Non-relativistic Hamiltonian operator 
with the coulombic interparticle interaction is used. Trial 
wavefunction to be corrected is a Slater determinant 
taken from the work of Emmanouilidou et al. [l4T | : 

^ T = e-^ ri+ar2 \l-ar 2 )-e-^ r2+ari \l-ar 1 ) , (7) 

where r\ and r 2 are the electron nucleus distances for the 
two electrons of the system and the numerical value of 



the parameter a is chosen as 0.65. This trial function 
does not satisfy the true nodal surface of the ls2s 3 S 
state of the Helium atom which is well known as r\ — 
Ti surface. Plots of the trial wavefunction at various 
cross sections show that the nodal surface of the trial 
wavefunction extends out of the true nodal surface. 

DMC time step is chosen as 0.0005 atomic time units 
and the data is collected for 80000 time steps. Initial 
number of plus and minus signed walkers and the number 
of points for the vacuum branchings are set to 500. The 
normalization of the trial wavefunction is set to 2600 and 
the number of walkers of each sign is stabilized around 
520, giving the r D value as 5.0. 

The non-relativistic true en erg y expectation value for 
this state is -2.1752 Hartrees [l5| and the energy expec- 
tation value of the used trial wavefunction is -2.1548 
Hartrees. Current correction scheme calculation gives 
the result as -2.1767(63) Hartrees for which the true 
value is in the statistical error interval. Computation 
time of the correction scheme calculation is 4.2 times 
shorter than the computation time to achieve the same 
precision with the comparison case calculation without 
any trial wavefunctions. Fixed-node DMC calculation 
using the nodes of the used trial wavefunction gives 
the energy value as -2.1626(8) Hartrees confirming the 
deviation of the trial wavefunction nodal surface from 
the true nodal surface r\ = r%- 



IV. DISCUSSION 

Application of the wavefunction correction scheme to 
the non fixed-node DMC increases the efficiency of the 
method significantly. Benchmark computations on the 
harmonic fermions and the helium atom show that the 
computation time for calculating the result within a cer- 
tain statistical error value decreases several times when 
the correction scheme is used. Improvement of the ef- 
ficiency depends on the ratio of the trial wavefunction 
normalization to the number of walkers and this param- 
eter of the method can be increased when the trial wave- 
function gets closer to the true fermionic wavefunction. 
Therefore, quality of the trial wavefunction is an impor- 
tant factor affecting the efficiency improvement observed 
in the correction scheme calculations. The mentioned 
ratio parameter value should be set to the highest value 
allowing a full correction which can be guessed consider- 
ing the quality of the trial wavefunction. This parameter 
value can be optimized using the fact that a deviation in 
the calculated expectation value starts to occur beyond 
the optimum value of the parameter. 

The nodal surfaces of the used trial wavefunctions de- 
viate from the true surfaces for the all cases except the 
harmonic fermions calculation in ID for which a permu- 
tation cell preserving node distortion is not possible. The 
correction scheme calculations yield the true energy ex- 
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pectation value despite the wrong nodal surfaces which 
proves the applicability of the wavefunction correction 
scheme for the non fixed-node QMC calculations. 

Fixed-node DMC is applicable to large systems since it 
eliminates the minus signed walkers by constraining the 
computation in a nodal region. This advantage of the 
fixed-node approximation disappears with the usage of 
the correction scheme because of the arising minus signed 
walkers as a result of the vacuum branchings. However, 
the wavefunction correction technique is suitable for non 
fixed-node DMC computations since the minus signed 
walkers are already needed for plus-minus cancellation 
methods. Correction scheme improves the large scale ap- 
plicability for such calculations as opposed to the case for 
the fixed-node calculations. 

Benchmark computations are carried out without us- 
ing the importance sampling transformation in order to 
better observe the sole effect of the wavefunction correc- 
tion technique. Importance sampling may be facilitated 
in the correction scheme calculations as described in the 
references [1, However, the guiding function should 
allow the walkers' passage from the boundaries of the 
chosen permutation cell for the application of the cur- 
rent boundary conditions when the fixed-node constraint 



is relaxed. A slightly modified form of the corrected trial 
wavefunction which does not vanish on the nodal sur- 
face of the original function may be used as the guiding 
function. 

A generally applicable method beyond the fixed-node 
approximation is very desirable for high accuracy elec- 
tronic structure calculations of relatively larger systems. 
Current non fixed-node QMC methods have exponen- 
tial scaling computation cost with increasing number of 
fcrmions and therefore not applicable to large systems. 
The fermion sign problem may not have a polynomial 
time solution with the classical computation techniques. 
However, fixed-node constraint may be relaxed by some 
sort of other approximate manners and the application of 
the non fixed-node computations may be widened by the 
usage of improvements like the wavefunction correction 
technique used in the current study. 
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